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Abstract. We discuss the fundamentals of the implicit moment method for 
Particle In Cell (PIC) simulation as presently implemented in the CELESTE3D 
code. We present the method in its fully electromagnetic and fully kinetic ver- 
sion. The application of the method is to problems with multiple temporal and 
spatial scales, common in all space, astrophysical and laboratory plasmas. 



1. Introduction 

The implicit particle in cell (PIC) method was developed as a general plasma 
simulation tool based on a fully kinetic approach (Masonlll981 : Denavit 198ll ) 



that did not rely on any physical approximation but only on the use of advanced 
numerical methods to reduce the cost of large scale kinetic simulations. 

In its simplest form the PIC method uses explicit time discretization meth- 
ods ( Birdsall and Langdon|[2004l ) . The equations of motion need the fields acting 



on the computational particles and the field equations need the moments com- 
puted from the particles. The great majority of PIC codes in use today address 
this problem relying on a explicit method. 

In the explicit method, the field equations are discretized with one in a num- 
ber of different explicit metho ds available (see iBirdsall and LangdonI 1^004 ) : 



Hocknev and Eastwood! ( 19881 ) for a review) and the equations of motion are 



usually discretized with the leap-frog algorithm ( Birdsall and LangdonI |2004 



Hocknev and Eastwood! [l988l ). The key point is that in the explicit method, 



the field equations need only the sources from the previous time cycle and the 
equations of motions need only the fields from the previous time cycle. Even 
though the equations remain coupled, no iteration is needed and the cycle be- 
comes a simple marching order where each block is applied after its predecessor 
and needs only information already available. This choice makes the explicit 
PIC very simple. As can be imagined this simplicity comes at a price. The 
explicit PIC approach is subject to three very restrictive stability constraints. 

First, the explicit discretization of the field equations requires that a Courant 
condition must be satisfied on the speed of light: 

cAt < Ax (1) 

Second, the explicit discretization of the equations of motion introduces 
a constraint related to the fastest electron response time, the electron plasma 
frequency: 

iOpeAt < Ax (2) 
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Note that, instead, thanks to the Boris algorithm for the motion in a magnetic 



field, the gyromotion introduces no stability constraints (jBirdsall and Langdon 
20041 ). 



Finally, the interpolation between grid and particles causes a loss of infor- 
mation and an aliasing instability called finite grid instability that results in an 
additional stability constraint: 

Ax < ?ADe (3) 

that requires the grid spacing to be of the order of the Debye length or smaller, 
the proportionality constant ^ being of order one and dependent on th e details 
of the scheme used (iB irdsall and LangdonI [20041 : iHocknev and Eastwoo d 1988). 

Thanks to ever faster computers, the explicit PIC method has been able to 
achieve remarkable results. But there are problems where multiple scales still 
prevent its use. When the interest is on large scales and slow processes other 
approaches are needed. The implicit PIC approach provides a viable alternative 
in such cases. 



2. Implicit PIC 



The implicit PIC method has been dev eloped along two dif ferent lines of in- 
vestigation: the direct implicit method ( Langdon et al. 19831 ) and the implicit 
moment method ( Brackbill and Forslu nd 19821) . Here we consider the implicit 
moment method. In the last few years, we have published a number of papers on 
the application o f CELESTE3D and a full description of the methods used has 
been provided bv I Lapenta et al" ( 20061 ). In the present work, we summarize the 
latest status of the implicit moment method as implemented in CELESTE3D. 

We consider first the time discretization. The solution is advanced in time 
with discrete steps. At, from the initial time, t*^ = to the final time = T. 
The generic quantity ^' at time step n {t = t"), is denoted with 

The equations of motion are discretized as ( Vu and Brackbill 19921 ): 



V -I- 



qsAt 



E^+^(x^+l/2) + ^n+1/2 ^ B^(x^+l/2; 



(4) 



where all quantities evaluated at intermediate levels are computed as: ^'"+^ = 
'^^{l—6)+'^"'~^^6. Note that the velocity equation is more conveniently rewritten 
as: 



,n+l/2 



n+9( n+l/2\ 



(5) 



where (3s = QpAt/rUp (independent of the particle weight and unique to a given 
species). For convenience, we have introduced hatted quantities obtained by 
explicit transformation of quantities known from the previous computational 
cycle: 



(6) 
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The transformation tensor operators a" are defined as: 

and represent a scaling and rotation of the velocity vector. 

The semi-discrete (continuous in space) temporal discretization to Maxwell's 
equations is written as: 

V X E"+^ + = 

c At 

c At c (8) 

V • E"+^ = 47r/)"+^ 

V • B" = V • B"+i = 0, 

The parameter 6 € [1/2, 1] is chosen in order to adjust the numerical dis- 
persion relation for electromagnetic waves (fo r 6 < 1/2, the algorithm is shown 
to be unstable ( Brackbill and Forslund|[l98^ )). We note that for 9 = 1/2 the 



scheme is second-order accurate in At; for 1/2 < < 1 the scheme is first-order 
accurate. 

The sources in Maxwell's equations necessitate information from the parti- 
cles: 

ps = Y.p=i qpW{yi - xp) 

(9) 

where the species is labelled by s and the sums are carried over all particles of 
a species Ng. The coupling with the particle equations of motion is evident. 

The fundamental problem to address in developing an implicit PIC method 
is the coupling between the equations of motion and the field equations for 
the presence of the time advanced electric field (but not magnetic field, that is 
used from the previous cycle, as no instability is introduced) in the equations 
of motion and for the appearance of the particle properties in the sources of 
the Maxwell equations. In both cases the coupling is implicit, so that the new 
particle properties need to be known before the fields can be computed and 
likewise the new fields need to be available before the new particle properties 
can be computed. 



3. Implicit Moment Method 

The implicit moment method removes the need for iterative methods and pro- 
vides a direct method to compute the advanced fields without first having to 
move the particles. The implicit moment method reduces the number of equa- 
tions that must be solved self-consistently to a set of coupled fluid moment and 
field equations. The solution of these equations implicitly, and the subsequent 
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solution of the particle equations of motion in the resulting fields, is stable and 
accurate. 

The coupling due to the implicit discretization of both field and particle 
equations is approximated, representing the sources of the field equations using 
the moment equations instead of the particle equations directly. Once the field 
equations are solved within this approximation, the rest of the steps can be 
completed directly without iterations: with the new fields, the particle equations 
of motion can be solved and the new current and density can be computed for 
the next computational cycle. 

The implicit moment i netho d formulation used he re is described in details 
bv lBrackbill and Forslundl (Il982l 1: iLaoenta et al.l (|2006l V The key step is to de- 



rive a suitable set of moment equations that can approximate the particle motion 
over a computational cycle. The approach followed in the present implementa- 
tion is based on a series expansion of the interpolation functions used to transfer 
information between grid and particles. The details of the siniple bu t demand- 
ing algebraic manipulations are provided bv IVu and Brackbilll ( 19921 ). the final 
answer being: 

where the following expressions were defined: 

J. = Ep qp%w{x - xp 
n. = Ep9pVpVpt^(x-xp 



(10) 



(11) 



with the obvious meaning, respectively, of current and pressure tensor based on 
the transformed hatted velocities. An effective dielectric tensor is defined to 
express the feedback of the electric field on the plasma current and density: 

..n QsPs ^^'^\ 
t^s = OLs (12) 

The expression (jlOp for the sources of the Maxwell's equations provide a 
direct and explicit closure of Maxwell's equations. When eq. (jlOp is inserted 
in eq. the Maxwell's equations can be solved without further coupling with 
the particle equations. This is the key property of the moment implicit method 
and allows the implicit moment PIC method to retain the once-through ap- 
proach typical of explicit methods and eliminates the need for expensive itera- 
tion procedures that would require to move the particles multiple times per each 
computational cycle. 



4. Stability 



The stability proper ties of the method described above have been studied ex- 
tensively in the past ( Brackbill and Forslund 19821 ). All the stability constraints 
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discussed above for the explicit method are removed. The imphcit particle mover 
removes the need to resolve the electron plasma frequency, and the implicit for- 
mulation of the field equations removes the need to resolve the speed of light. 

The time step constraints are replaced by an accuracy limit arising from the 
derivation of the fluid moment equations using the series expansi on. This limit 



restricts the mean particle motion to one grid cell per time step (|Brackbill and Forslund 
19821 ^. i.e. 



Vth,e^t/Ax < 1, (13) 

Th e finite grid instability limit f or the explicit method, Ax < <rA/5e is replaced 
by (jBrackbin and Forslundlll982l 1 



Ax/ At < qvth,e, (14) 

that allows large grid spacings to be used when large time steps are taken. The 
gain afforded by the relaxation of the stability limits is two-fold. 

First, the time step can far exceed the explicit limit. In a typical plasma the 
electron plasma frequency is far smaller than the time scales of interest and its 
accurate resolution is not needed. Within the current approach, the processes 
developing at the sub-At scale are averaged and their energy is damped by a 
numerically-enhanced La ndau darnping . In other approaches, such as the gyroki- 
netic or hybrid approach (lLipatovll2002l ). such processes are completely removed 



and the energy channel towards them is interrupted, removing for example the 
possibility to exchange energy between sub-At fluctuations and particles. In 
the implicit approach, instead, the sub-At scales remain active and the energy 
channel remains open. This is a crucial feature to retain a full kinetic approach. 
Furthermore, when additional resolution of the smallest scales is needed, the im- 
plicit method can access the same accuracy of the explicit method simply using 
a smaller time step and grid spacing. This feature is not accessible to reduced 
models, e.g. gyroaveraged methods, that remove the small scales entirely. 

Second, the grid spacing can far exceed the Debye length. Often the scales 
of interest are much larger than the Debye length. The ability to retain a full 
kinetic treatment without the need to resolve the Debye length results in a much 
reduced cost for the implicit PIC method. 



5. Conclusions 

The implicit moment PIC method described above is implemented in the CE- 
LESTE3D code. The C ELESTE3D code was origina lly conceived for the nu- 
merical tokamak project ( Brackbill and Lapental [l994l ) but has found its main 



application in space physics. 

Four types of tests have been conducted to verify and validate CELESTE 
in full 3D cartesian geometry and in reduced geometries in 2D and ID: 1) well 
known benchmarks in cluding shocks, the We ibel instability, Landau damping 



and io n acoustic waves (|Vu and Brackb ill 1992); 2 ) the G EM challenge (jRicci et al 



2002bl lal) and the Newton challenge (iBirn et aD l2005l ^ : 3) study of reconnec- 



tion in systems with low betas, investigating the reconnection process both at 
the macroscopic and microscopic le vel obtaining agreement with the explicit 
PIC code NPIC (jRicci et al.l l2004al ) : 4) 3D stability study of a current sheet 
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Figure 1. 3D Simulation of the magnetotail requiring 1 day of CPU time 
on CELESTE3D or 800,000 years on a explicit code 



equilibrium, compared with satellite observations obtain ed from the CLUS- 
TER and GEOTA IL mission (.Lapenta and Brackbill ,2002i : Lapenta et al. 20031 : 
iRicci et alJl2004bl 1. 

As an example of this last case, we show here the results of an actual 
simulation co nducted in 3D of a system initially in the magnetotail equilibrium 
described by iBirnI ( 19871 ). The physics developing in the simulation includes 
first the growth and saturation of the lower hybrid drift instability, followed 
by the onset of reconnection and current flapping leading to the macroscopic 
restructuring of the topological configuration of the magnetotail. The physics 
steps of the process are described in de atils by iLapenta and Brackbilll ( 20021 . 
20001 1: iLaoenta et al.l ([2003 1 : iRicci et al.1 (|2004bl ^. 

The fundamental consideration of interest here is the efficacy of the simula- 
tion approach. In a typical explicit run, the time step would have to be selected 
according to the stability constraint of ujpeAt < 2. In a typical magnetotail 
case, the electron plasma frequency is of the order of ujpe ~ 5 • 10^s~^, and the 
ion plasma frequency is of the order of Upi ~ lO^s^^, both smaller than the 
smallest scale of interest for this problem which is the lower hybrid frequency 
range of wlh ~ 10^s~^. In a implicit simulation, we can select the time step to 
the ion plasma frequency, still resolving accurately the lower hybrid range, but 
saving two orders of magnitude compared with the explicit case that instead is 
needlessly resolving the electron plasma scale. 

A similar gain occurs in each spatial direction. In a explicit simulation, for 
stability the grid spacing needs to satisfy Ax/Xde < The Debye length in the 
magnetotail is of the order of 100m. But the smallest scales of interest in the 
present case are in the range between the electron (10km) and the ion (100km) 
inertial scales or gyroscales. In our simulation we set the resolution at 10km, 
saving two orders of magnitude in each spatial direction but still resolving the 
important scales. 

Counting all savings, we have saved 2 orders of magnitude in each spatial 
direction and in time (in total 8 orders of magnitude) , without loosing the deatils 
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of the scales of interest. The savings means that on the same computer, an 
implciti simulation can compute in one day what an explicit simulation would 
require nearly 800,000 years. 
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